Context: Atmospheric Heating¶
How is the solar atmosphere heated?
How is the energy transported from the top of the convection zone / photosphere upwards?
How is the Energy Deposited?
Problem: MHD Wave Generation¶
"The challenge remains to understand how, and where, these waves are generated and dissipated in the solar atmosphere, [...]" - (McIntosh et al. 2011)
The Photosphere as a Source for MHD Waves¶
The dynamic photosphere with embedded magnetic field provides many potential ways of driving MHD waves.
- 'Buffeting' from convective motions
- Global oscillations e.g. (Ulrich 1970; Leighton, Noyes, and Simon 1962).
- Spiralling in downdrafts e.g. (Wedemeyer-Böhm et al. 2012; Bonet et al. 2010; Bonet et al. 2008)
Numerical Simulations of the lower solar atmosphere¶
The Code¶
The code used is the Sheffield Advanced Code (SAC) (Shelyag, Fedun, and Erdélyi 2008).
SAC is based on the Versatile Advection Code (VAC) (Tóth 1996).
SAC simulates perturbations on a static background, using a CD4 solver with hyper-diffusion and hyper-viscosity terms to stabalise the solution.
This makes it well-suited to solving wave perturbations on top of a highly stratified background such as the solar atmosphere.
The Model¶
To simulate wave excitation in the photosphere a numerical model of the solar atmosphere is needed.
Hydrostatic background from the VAL III model C (Vernazza, Avrett, and Loeser 1981):
import pysac.mhs_atmosphere as atm
#Read in the VAL3C model
empirical_data = atm.hs_atmosphere.read_VAL3c_MTW(MTW_file=False)
# Create a Z array at the interpolated resolution and interpolate.
ZZ = u.Quantity(np.linspace(empirical_data['Z'][0], empirical_data['Z'][-1], 128), unit=empirical_data['Z'].unit)
table = atm.hs_atmosphere.interpolate_atmosphere(empirical_data, ZZ, s=0)
# Create a figure and make space for the axes.
fig, ax = plt.subplots(gridspec_kw={'right':0.77, 'left':0.16, 'bottom':0.13}, figsize=(13,8))
# Shortcut all the Mm conversion.
Z = empirical_data['Z'].to(u.Mm)
lrho, = ax.plot(Z, empirical_data['rho'].quantity.si, 'x', color='blue')
lrho_i, = ax.plot(ZZ.to(u.Mm), table['rho'].quantity.si, color='blue')
ax2 = ax.twinx()
lp, = ax2.plot(Z, empirical_data['p'].to(u.Pa), 'x', color='green')
lp_i, = ax2.plot(ZZ.to(u.Mm), table['p'].to(u.Pa), color='green')
ax3 = ax.twinx()
ax3.spines["right"].set_position(("axes", 1.2))
lt, = ax3.plot(Z, empirical_data['T'].to(u.K), 'x', color='red')
lt_i, = ax3.plot(ZZ.to(u.Mm), table['T'].to(u.K), color='red')
# Set primary axes properties and labels
ax.semilogy()
ax.set_ylabel(r"Density [{}]".format(lrho._yorig.unit._repr_latex_()))
ax.set_xlabel(r"Height [{}]".format(lrho._xorig.unit._repr_latex_()))
ax.set_xlim(Z[0].value-0.1, Z[-1].value+0.1)
# Pressure Axis
ax2.semilogy()
ax2.set_ylabel(r"Pressure [{}]".format(lp._yorig.unit._repr_latex_()))
# Temp axis
ax3.set_ylabel(r"Temperature [{}]".format(lt._yorig.unit._repr_latex_()))
ax.set_xlim([-0.02,1.62])
ax3.set_ylim([3500,7500])
# Set the colours for the ticks and the labels.
ax.tick_params(axis='y', colors=lrho.get_color())
ax2.tick_params(axis='y', colors=lp.get_color())
ax3.tick_params(axis='y', colors=lt.get_color())
ax.yaxis.label.set_color(lrho.get_color())
ax2.yaxis.label.set_color(lp.get_color())
ax3.yaxis.label.set_color(lt.get_color())
fig
The Magnetic Flux Tube¶
The magnetic field model mimics a Magnetic Bright Point (MBP) and follows (Gent, Fedun, Mumford and Erdélyi 2013) and is constructed as a self-similar (Schlüter and Temesváry 1958; Low 1980), cylindrically symmetric, expanding field.
slc = yt.SlicePlot(ds, 'x', 'density', origin='lower-center-domain', axes_unit='Mm')
slc.set_figure_size(9)
slc.set_cmap('all', 'viridis')
slc.set_font_size(20)
slc.set_zlim('all', 0, 2.5e-4)
seed_points = np.zeros([11,2]) + 1.52
seed_points[:,0] = np.linspace(-0.99, 0.95, seed_points.shape[0], endpoint=True)
slc.annotate_streamlines('mag_field_y', 'mag_field_z', field_color='magnetic_field_strength',
plot_args={'start_points':seed_points, 'density':15, 'cmap':'Blues', 'linewidth':2,
'norm': matplotlib.colors.LogNorm(*ds.all_data().quantities.extrema("magnetic_field_strength"))})
#force render
slc.save('/tmp/test.png')
#use the raw Figure for transparent bg
slc.plots['density'].figure
Driving Waves in the Simulation Domain¶
$$ \vec{V}(x,y,z,t) = \vec{F}(x,y,z) \ e^{-\left(\frac{x^2}{\Delta x^2} + \frac{y^2}{\Delta y^2} + \frac{z^2}{\Delta z^2}\right)} \sin \left(2\pi \frac{t}{P}\right) $$
Identifying Waves from Broadband Drivers¶
Photospheric drivers excite multiple wave modes simulatenously.
How to quantify the relative strengths of the different modes from different drivers?
Assume uniform media:
- Locally decompose perturbations into Fast, Slow and Alfvén modes.
- Compare the percentage wave energy in each mode.
Relationship between Mode and Velocity Perturbation¶
In a high-$\beta$ plasma, and assuming parallel propagation ($k_\parallel >> k_\perp$) the three MHD wave modes of a uniform plasma pertub different components of velocity:
$$ \frac{|v_\parallel|}{|v_\perp|+|v_\parallel|} = \frac{1}{\frac{k_\parallel}{k_\perp}+1}\\ \ \\ |v_\parallel| >> |v_\perp| $$
$$ |v| = |v_\phi| $$
$$ \frac{|v_\parallel|}{|v_\perp|+|v_\parallel|} = \frac{1}{\frac{k_\perp}{k_\parallel}+1}\\ \ \\ |v_\parallel| << |v_\perp| $$
Decomposing Velocity Pertubations¶
xc = domain['xmax']/2
yc = domain['ymax']/2
ti = 0
n = 100
surf_seeds = []
for theta in np.linspace(0, 2 * np.pi, n, endpoint=False):
surf_seeds.append([tube_r * np.cos(theta + 0.5 * ti) + xc,
tube_r * np.sin(theta + 0.5 * ti) + yc, domain['zmax']])
seeds = np.array(surf_seeds)
#Add axes:
axes, outline = mpf.add_axes(np.array(zip(ds.domain_left_edge,ds.domain_right_edge)).flatten()/1e8)
#Add seed points to plot:
seed_points = mlab.points3d(seeds[:,0], seeds[:,1], seeds[:,2],
color=(0.231, 0.298, 0.752), scale_mode='none',
scale_factor=1.5)
mlab_view(fig.scene)
field_lines = SeedStreamline(seed_points = np.array(seeds))
bmag.add_child(field_lines)
field_lines.actor.mapper.scalar_visibility = False
field_lines.actor.property.color = (0,0,0)
field_lines.actor.property.line_width = 1.5
mlab_view(fig.scene)
pd_seeds = ttf.make_circle_seeds(100, 60, **domain)
fieldlines, surface = ttf.create_flux_surface(bfield.outputs[0], pd_seeds)
surface.output.lines = None
flux_surface = mlab.pipeline.surface(surface.output)
flux_surface.actor.mapper.scalar_visibility = False
flux_surface.actor.property.color = (0.8,0.8,0.8)
#flux_surface.actor.property.line_width = 0
mlab_view(fig.scene)
axes.visible = False
outline.visible = False
flux_surface.actor.property.edge_visibility = True
mlab_view(fig.scene, azimuth = 90, elevation = 75, distance=80, focalpoint=[63, 120, 110], aa=20)
poly_norms = ttf.make_poly_norms(surface.output)
normvec = mlab.pipeline.glyph(poly_norms.output)
normvec.glyph.glyph_source.glyph_source = normvec.glyph.glyph_source.glyph_dict['arrow_source']
normvec.glyph.glyph.scale_mode = 'data_scaling_off'
normvec.glyph.glyph.color_mode = 'color_by_scale'
normvec.glyph.glyph.scale_factor = 5
normvec.glyph.glyph_source.glyph_position = 'tail'
mlab_view(fig.scene, azimuth=85, elevation=80, distance=50, focalpoint=[63, 120, 110], aa=20)
MHD Waves excited by Different Photospheric Drivers¶
Mumford, S. J., Fedun, V., Erdélyi, R. - The Astrophysical Journal - January 2015 - Volume 799, Issue 1
Generation of Magnetohydrodynamic Waves in Low Solar
Atmospheric Flux Tubes by Photospheric Motions
Driving Waves in the Simulation Domain¶
$$ \vec{V}(x,y,z,t) = \vec{F}(x,y,z) \ e^{-\left(\frac{x^2}{\Delta x^2} + \frac{y^2}{\Delta y^2} + \frac{z^2}{\Delta z^2}\right)} \sin \left(2\pi \frac{t}{P}\right) $$
(Wedemeyer-Böhm and Rouppe van der Voort 2009)
Archmedian Spiral¶
$$ F_x = A \frac{B_Ax}{x^2 + y^2} \frac{y}{\sqrt{x^2 + y^2}},\\ F_y = - A \frac{B_Ay}{x^2 + y^2} \frac{x}{\sqrt{x^2 + y^2}}. $$
Logarithmic Spiral¶
$$ F_x = A \frac{\cos(\theta + \phi)}{\sqrt{x^2 + y^2}},\\ F_y = - A \frac{\sin(\theta + \phi)}{\sqrt{x^2 + y^2}},\\ $$ where, $\theta = \tan^{-1}\left(\frac{y}{x}\right),\ \phi = \tan^{-1}\left(\frac{1}{B_L}\right)$
Analysis and Results¶
Wave Flux¶
Calculate wave energy flux from (Leroy 1985).
$$
\vec{F}_{\text{wave}} \equiv \widetilde{p}_k \vec{v} + \frac{1}{\mu_0} \left(\vec{B}_b \cdot \vec{\widetilde{B}}\right) \vec{v} - \frac{1}{\mu_0}\left(\vec{v} \cdot \vec{\widetilde{B}} \right) \vec{B}_b
$$
Decompose onto the flux surfaces in the same manner as velocity.
Time-Distance Diagrams¶
Vertical Driver¶
figures['vert']
(Mumford, Fedun, and Erdélyi 2015)
Horizontal Driver¶
figures['horiz']
(Mumford, Fedun, and Erdélyi 2015)
Logarithmic Spiral¶
figures['Slog']
(Mumford, Fedun, and Erdélyi 2015)
Conclusion 1: Driver Comparison¶
from flux_comparison import make_flux_bar_chart, get_averages
Favgs = get_averages(ch4.data_dir)
fig = make_flux_bar_chart((12,10), ch4.data_dir)
fig
(Mumford, Fedun, and Erdélyi 2015)
Effects of Logarithmic Spiral Expansion Factor¶
Mumford, S. J. and Erdélyi, R. - Monthly Noticies of the Royal Astronomical Society - March 2015 - Volume 449 Issue 2.
Photospheric Logarithmic Velocity Spirals as MHD Wave Generation Mechanisms
Expansion Factor¶
$$ B_L^{-1} = 6.4 \pm 1.6 \equiv B_L = 0.15 $$
More Time-Distance Diagrams¶
Spiral Expansion Factor $B_L=0.015$¶
bl_figures['B0015']
(Mumford and Erdélyi 2015)
Spiral Expansion Factor $B_L=0.15$¶
bl_figures['B015']
(Mumford and Erdélyi 2015)
Spiral Expansion Factor $B_L=1.5$¶
bl_figures['B15']
(Mumford and Erdélyi 2015)
Conclusion 2: Spiral Expansion Factor Comparison¶
import os
int_periods = np.floor(600./cfg.period)*180
def calc_avgs(tube_r):
AvgsP = np.zeros([3,len(BL)])
for i, bl in enumerate(BL):
cfg.exp_fac = bl
times = np.load(os.path.join(cfg.data_dir, 'Times_{}.npy'.format(cfg.get_identifier())))
max_index = np.argmin(np.abs(int_periods - times))
Fpar, Fperp, Fphi = map(np.load, ph.glob_files(cfg, tube_r, 'LineFlux*Fp*npy'))
#Fpar, Fperp, Fphi = map(np.load, ph.glob_files(cfg, tube_r, '*vp*npy'))
Fpar[np.abs(Fpar)<1e-5], Fperp[np.abs(Fperp)<1e-5], Fphi[np.abs(Fphi)<1e-5] = 0., 0., 0.
Fpar, Fperp, Fphi = Fpar[:max_index,:], Fperp[:max_index,:], Fphi[:max_index,:]
Ftot2 = (Fpar**2 + Fperp**2 + Fphi**2)
Fpar2, Fperp2, Fphi2 = np.array([Fpar, Fperp, Fphi])**2
FparP, FperpP, FphiP = (Fpar2/Ftot2)*100, (Fperp2/Ftot2)*100, (Fphi2/Ftot2)*100
FparP = np.mean(np.ma.masked_array(FparP, np.isnan(FparP)))
FperpP = np.mean(np.ma.masked_array(FperpP, np.isnan(FperpP)))
FphiP = np.mean(np.ma.masked_array(FphiP, np.isnan(FphiP)))
AvgsP[:, i] = FparP, FperpP, FphiP
return AvgsP
fig, axs = plt.subplots(nrows=3, figsize=(11,11), sharex=True)
titles = ["Flux Surface Radius $=156$ km", "Flux Surface Radius $=468$ km", "Flux Surface Radius $=936$ km"]
tubes = []
for i, ax in enumerate(axs):
AvgsP = calc_avgs(cfg.tube_radii[i])
tubes.append(AvgsP)
ax.semilogx()
ax.plot(BL, AvgsP[0], 'o', label=r"$F_\parallel^2$", mew=0, ms=10)
ax.plot(BL, AvgsP[1], '_', label=r"$F_\perp^2$", mew=2, ms=10)
ax.plot(BL, AvgsP[2], 'x', label=r"$F_\phi^2$", mew=2, ms=10)
ax.set_ylabel("% Square Wave Flux")
ax.set_title(titles[i])
ax.xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
ax.xaxis.set_ticks(BL)
ax.set_xlim([0.01, 2.01])
ax.set_ylim([0, 85])
err = np.array([-1*(0.15-1/(6.4-1.6)), 0.15+1/(6.4+1.6)])
ax.fill_betweenx(np.linspace(-5,105), err[0], err[1], alpha=0.3, color='green', linewidth=0)
axs[0].legend(bbox_to_anchor=(1.20, 1.05))
axs[-1].set_xlabel("Logarithmic Spiral Expansion Factor ($B_L$)")
plt.tight_layout()
fig
(Mumford and Erdélyi 2015)
Maintaining Consistent Energy Input¶
The amplitude is calculated using:
$$
A^2 \propto \frac{1}{P}
$$
with $A=10$ ms$^{-1}$ and $P=240$ s as the reference point.
| Period [seconds] | Amplitude [ms$^{-1}$] |
| $30.0$ | $20\sqrt{2}$ |
| $60.0$ | $20$ |
| $90.0$ | $20\sqrt{\frac{2}{3}}$ |
| $120.0$ | $10\sqrt{2}$ |
| $150.0$ | $4\sqrt{10}$ |
| $180.0$ | $\frac{20}{\sqrt{3}}$ |
| $210.0$ | $20\sqrt{\frac{2}{7}}$ |
| $240.0$ | $10$ |
| $270.0$ | $\frac{20}{3}\sqrt{2}$ |
| $300.0$ | $4\sqrt{5}$ |
Even More Time-Distance Diagrams¶
$P=30$ $s$ $A=20\sqrt{2}$ ms$^{-1}$¶
p_figures['30.0']
$P=150$ $s$ $A=4\sqrt{10}$ ms$^{-1}$¶
p_figures['150.0']
$P=300$ $s$ $A=4\sqrt{5}$ ms$^{-1}$¶
p_figures['300.0']
Conclusion 3: Period Comparison¶
from period_amps import periods, sim_params
sim_params = sim_params[:10]
periods = np.array(periods[:10])
fig, axs = plt.subplots(nrows=3, figsize=(11,11), sharex=True)
titles = ["Flux Surface Radius $=156$ km", "Flux Surface Radius $=468$ km", "Flux Surface Radius $=936$ km"]
tubes = []
for i, ax in enumerate(axs):
AvgsP = ph.get_all_avgs(cfg, cfg.tube_radii[i], sim_params)
tubes.append(AvgsP)
ax.plot(periods, AvgsP[0], 'o', label=r"$F_\parallel^2$", mew=0, ms=12)
ax.plot(periods, AvgsP[1], '_', label=r"$F_\perp^2$", mew=2, ms=12)
ax.plot(periods, AvgsP[2], 'x', label=r"$F_\phi^2$", mew=2, ms=12)
ax.set_ylabel("% Square Wave Flux")
ax.set_title(titles[i])
ax.xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
ax.xaxis.set_ticks(periods)
ax.set_xticklabels(["{}\n[{:n}]".format(p, 600//p) for p in periods])
ax.set_ylim([10, 75])
ax.set_xlim([25, 305])
axs[-1].set_xlabel("Period [s] \n [Number of periods averaged]")
#axs[0].legend(loc=2)
axs[0].legend(bbox_to_anchor=(1.21, 1.05))
plt.tight_layout(h_pad=0.1)
fig
Conclusions¶
- Different spiral driver profiles do not have a large impact on the excited waves.
- Spirals with small expansion factors only excite $40-60$% Alfvén mode.
- Spirals with large expansion factors exicte mostly fast mode.
- Period only changes excited flux by $10-20$%
My Publications¶
Mumford, S. J. and Erdélyi, R. - Monthly Noticies of the Royal Astronomical Society - March 2015 - Volume 449 Issue 2.
Photospheric Logarithmic Velocity Spirals as MHD Wave Generation Mechanisms
Mumford, S. J., Fedun, V., Erdélyi, R. - The Astrophysical Journal - January 2015 - Volume 799, Issue 1
Generation of Magnetohydrodynamic Waves in Low Solar
Atmospheric Flux Tubes by Photospheric Motions
The SunPy Community, Mumford, S. J., Christe, S., Pérez-Suárez, D., et. al - Computational Science and Discovery - January 2015 - Volume 8 Issue 1.
SunPy: Python for Solar Physics
Freij N., Scullion E. M., Nelson C. J., Mumford S. J., Wedemeyer S., and Erdélyi R. - The Astrophysical Journal - July 2014 - Volume 791, Issue 1, p.61
The Detection of Upwardly Propagating Waves Channeling Energy from the Chromosphere to the Low Corona
Gent, F. A., Fedun, V., Mumford, S. J., Erdélyi, R. - Monthly Notices of the Royal Astronomical Society - October 2013 - Volume 435, Issue 1, p.689-697
Magnetohydrostatic equilibrium - I. Three-dimensional open magnetic flux tube in the stratified solar atmosphere
Nelson, C. J., Doyle, J. G., Erdélyi, R., Huang, Z., Madjarska, M. S., Mathioudakis, M., Mumford, S. J., Reardon, K - Solar Physics - April 2013 - Volume 283, Issue 2, p.307-323.
Statistical Analysis of Small Ellerman Bomb Events